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We investigate the tricritical scaling behavior of the two-dimensional spin-1 Blume-Capel model by using 
the Wang-Landau method measuring of the joint density of states for lattice sizes up to 48 x 48 sites. We find 
that the specific heat deep in the first-order area of the phase diagram exhibits a double-peak structure of the 
Schottky-like anomaly appearing with the transition peak. The first-order transition curve is systematically de¬ 
termined by employing the method of field mixing in conjunction with finite-size scaling, showing a significant 
deviation from the previous data points. At the tricritical point, we characterize the tricritical exponents through 
finite-size-scaling analysis including the phenomenological finite-size scaling with thermodynamic variables. 
Our estimation of the tricritical eigenvalue exponents, yt = 1.804(5), y g = 0.80(1), and yh = 1.925(3), 
provides the first Wang-Landau verification of the conjectured exact values, demonstrating the effectiveness of 
the density-of-states-based approach in finite-size scaling study of multicritical phenomena. 
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I. INTRODUCTION 

The Wang-Landau (WL) sampling method mm directly 
estimates the density of states through random walk in energy 
space. Because of its capability of dealing with complex en¬ 
ergy landscape together with the flexibility for applications, it 
has been widely used in different areas of physics and chem¬ 
istry, including protein folding J3] |4), fluid simulations 0, 
random spin systems ED, and also quantum systems cziii. 
Particularly for study of phase transitions, the WL method 
suggests an efficient way to overcome the issue of slow dy¬ 
namics in the conventional Monte Carlo simulations. Reduc¬ 
ing tunneling time and critical slowing down in the first- and 
second-order transitions has been a long-standing subject in 
the advances of the Monte Carlo methods which include, for 
instance, the cluster algorithms mm, multicanonical en¬ 
semble mum parallel tempering mm, and histogram 
reweighting technique ED- In the WL method, with the den¬ 
sity of states being accurately estimated, one can immediately 
access thermodynamic quantities at any temperatures across 
phase diagram, indicating its potential for study of critical 
phenomena (for instance, see liT6ll22l ). In this paper, we focus 
on the tricritical phenomena and examine the effectiveness of 
the Wang-Landau method in the finite-size-scaling analysis of 
the tricritical behavior in two dimensions. 

A tricritical point at which the nature of phase transition 
changes from first order to second order has been observed in 
a variety of systems (23ll : for instance, multicomponent flu¬ 
ids, metamagnets, i He- 4 He mixtures (24j, and also recently 
ultracold quantum gases f25j . Interestingly, the upper critical 
dimension for the Ising tricritical behavior is lowered to three, 
and thus in two dimensions, the tricritical scaling exponents 
become different from the classical ones 1261127]. The tricrit¬ 
ical universality in two dimensions was first conjectured from 
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the dilute Potts model (28l[30l and established by the confor¬ 
mal invariance argument ED- On the other hand, large efforts 
with advanced numerical methods on different models have 
been also devoted to precisely calculate the tricritical eigen¬ 
value exponents, namely the thermal exponent y t , the next- 
to-leading thermal exponent y g , and the magnetic exponent 
yh- The Monte Carlo renormalization group (MCRG) calcu¬ 
lation was performed for the Ising antiferromagnet (AFM) and 
the Blume-Capel (BC) model (32]; the transfer-matrix method 
was applied to the BC model (33ll : the metropolis algorithm 
with the histogram reweighting (HR) technique was used for 
the spin fluid and the BC model l34l . While the WL method 
was first applied to the BC model in Ref. im, the tricritical 
scaling exponents still remain unexplored in the same method. 
The tricritical eigenvalue exponents estimated from these pre¬ 
vious calculations are listed in Table U 

Here, by using the Wang-Landau method, we approach 
the tricriticality of the two-dimensional spin-1 Blume-Capel 
model from the side of the first-order phase transitions. The 
joint density of states is measured for systems with sizes up 
to 48 x 48 sites, allowing an accurate picture of the first-order 
transitions and tricritical scaling behavior. First, at a large 
crystal field, we find out a double-peak structure of the spe¬ 
cific heat where the Schottky-like anomaly appears together 
with the first-order transition peak. It turns out that our large- 
scale calculations are crucial to reveal this anomalous struc¬ 
ture. Second, we systematically determine the first-order tran¬ 
sition curve which provides significant deviations from the 
implicit line of the few previous data points. Finally, we char¬ 
acterize the tricritical exponents through finite-size-scaling 
analysis within the Wang-Landau framework. The excellent 
data-curve collapses in the finite-size scaling accurately deter¬ 
mine the three tricritical eigenvalue exponents, providing the 
first Wang-Landau verification of the conjectured exact values 
of the tricritical exponents in two dimensions. 

This paper is organized as follows. Section [IT] defines the 
Hamiltonian of the BC model and provides the numerical de¬ 
tails of our simulations. In Sec. [II] we also briefly describe 
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Numerical method 

Model 

Vt 

Vg 

Vh 

MCRGf32l 

Ising AFM, BC 

1.80(2) 

0.84(5) 

1.93(1) 

Transfer Matrix l33l 

BC 

1.75(3) 

0.80(1) 

1.90(5) 

Metropolis+HR (34l 

Spin fluid, BC 

1.80(1) 

0.83(5) 

1.93(1) 

WL (this work) 

BC 

1.804(5) 

0.80(1) 

1.925(3) 

The exact conjectures 12814311 

9/5 

4/5 

77/40 


TABLE I. Numerical estimations of the tricritical eigenvalue expo¬ 
nents in two dimensions. The conjectured exact exponents are given 
for comparison. 


the mixing-field method with which we proceed for the char¬ 
acterization of the tricritical behavior of the BC model. In 
Sec. Ill we present the specific-heat anomaly observed with 
the first-order transitions, and we provide the determination 
of the first-order-transition line with the details of the mixing- 
field analysis ||34| - (36l that we employ to examine phase co¬ 
existence and to locate the tricritical point. In Sec. IV we 


provide the finite-size-scaling analysis to estimate the tricrit¬ 
ical exponents, including the one with the probability distri¬ 
bution associated with the field mixing and the phenomeno¬ 
logical finite-size scaling with thermodynamic quantities. In 
Sec.[V] summary and outlooks are given. 


II. MODEL AND METHODS 

A. Grand canonical formulation of the Blume-Capel model 

The spin-1 Blume-Capel model in square lattices that we 
consider can be described by the Hamiltonian 

H = s i s j + - h^2 s i, ( 1 ) 

<».j> » * 


where spin s* at site i can take a value of +1, —1, or 0, and 
J and A denote the ferromagnetic coupling and crystal field 
causing spin anisotropy, respectively. The summation £ 
runs over all pairs of nearest-neighbor spins. The coupling J 
is set to be unity to define unit energy scale. We only consider 
the case of zero external magnetic field, h = 0, in the calcula¬ 
tions. The system size is denoted by L representing L d lattice 
sites where the dimension d = 2 for our square lattices. For 
the numerical implementation, we write the partition function 
in a grand canonical form as 

Zl(P, = r ( S > N > N e MPE), ( 2 ) 

E,N 

where f3 denotes the inverse temperature l/fc^T, and the fu- 
gacity z is given as z = exp (—fi) with p = (3 A. The Boltz¬ 
mann constant kg is set to unity for simplicity. The variables 
E = /( . : p. SiSj and N = sf represent the kinetic energy 

and number of nonzero spins, respectively. The joint density 
of states F(E, N) is to be given by the WL sampling. 


B. Direction of scaling fields 

We employ the method of field mixing ll34U36l to describe 
the asymmetry of phase transition undergoing in the Blume- 
Capel model and the scale invariance at the tricritical point. 
The formulation in Eqs. G and (|2]i suggests the temperature 
T, crystal field A, and magnetic field h as a natural choice 
of fields to describe the phase diagram. While the scaling 
direction associated with h is orthogonal to the T — A (or 
j3 — n) plane because of the Ising symmetry, there is no such 
symmetry for the other two fields. Thus, for instance in order 
to study the tricritical behavior, one may consider the linear 
combinations of /3 and // to describe the relevant scaling fields 
as 


A = {fj, - m) + r(/3 - 0 t ), (3) 

9 = (P - Pt) + s{p- nt), (4) 

where p t and p t are the values at the tricritical point, and r and 
s are the mixing parameters. The scaling field g is tangent to 
the coexistent curve while the direction of A is not restricted. 
Accordingly, one can also write down the two relevant vari¬ 
ables conjugate to the scaling fields as 

Q =— 1 —(n — se), (5) 

1 — rs 

£ = —( e ~ rn )’ ( 6 ) 
1 — rs 

where n = L~ d N and e = L~ d E, satisfying the requirement 
( X) = L~ d d In Zg/dx for scaling field x. Note that the mix¬ 
ing parameters r and s are system-specific quantities, and thus 
the scaling fields and their corresponding conjugate variables 
can exhibit system-size dependence. 

In the vicinity of the tricritical point, the finite-size-scaling 
ansatz for the limiting probability distribution function of the 
scaling fields and their conjugate variables is written as 

P L a p L (a^ 1 L d ~ Vt Q, a^L^S, a^L^m, 

a t L Vt X, a g L Va g, a h L Vh h), (7) 

where p/\ is a universal scaling function, and a’s are nonuni- 
versal factors (for more details, see Refs. iEHSI)- Precisely 
at the tricritical point, the probability distribution function be¬ 
comes 


P L oc Pl^ 1 L d ~ yt Q, a~ 1 L d ~ Va £,a^ 1 L d ~ yh m), (8) 

where p} is universal and scale invariant, which allows use to 
measure the tricritical exponents y ’s from the finite-size scal¬ 
ing for systems with different sizes. The probability distribu¬ 
tion function of the field-conjugate variables can be estimated 
from the histogram accumulating the occurrence of (E, Ap¬ 
points in the discrete bins of Q and £ with weighting factor 
T{E,N)z n exp{/3E). 


C. Numerical aspects of the Wang-Landau sampling 

Our numerical estimation of I N) follows the standard 
WL algorithm EH except that our random walk needs to be 
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performed in the two dimensional parameter space of E and 
N. Initially, the density of states Y(E. N) is set to be unity, 
and the random walk proceeds by trying out a new random 
value for a spin randomly chosen in the lattices. The new trial 
spin would move the energy from {E\, Ni) to (E?, N%) in the 
two dimensional parameter space, and then this spin update is 
accepted with the transition probability 

P [(£„ JV.) -» (E 2 , <V 2 )] = min l) (9) 

for the sake of the importance sampling of T. In ev¬ 
ery trial of spin updates, the current energy state (E, N) 
is recorded in T and the histogram H of visited states as 
In T(E, N ) —> In T(E, N ) + In / with the modification factor 
/ and H(E,N) —> H(E,N) + 1, respectively. These WL 
procedures continues until the histogram becomes flat enough 
and then are restarted with a reduced modification factor and 
with resetting the histogram as H = 0. 

The modification factor / is initially given as In / = 1 and 
scaled down as when restarting. The flatness criterion 
for the histogram is set to be 95% for L = 16 and 90% for 
L = 20 and 24, and it is lowered to 80% for L > 32. To 
avoid accidental satisfaction of the flatness criterion, the num¬ 
ber of Monte Carlo steps (MCS) between successive flatness 
inquiries is set to be the same as the number of available en¬ 
ergy states of (E, N ), where unit MCS is defined as L d trials 
of a single spin update. The actual flatness inquiry interval is 
about 10 5 MCSs for L = 16 and increases to 10 7 for L = 48. 
The stopping criterion for the modification factor is given as 
In / < 10 -8 for L < 32, a less stringent 10~ 7 for L = 40, 
and 10" 6 for L = 48. 

The main difficulty encountered in these procedures comes 
from the increased dimensionality encoded in the importance 
sampling with spin update trial, causing very long computa¬ 
tion times. Compared to the usual case with a single energy 
parameter, many more spin update trials are needed to cover 
the two-dimensional space of (E, N) since one spin update 
trial can visit only one energy state. For L = 48, the size 
of (E, iV)-space is enlarged to be about 10', while the cor¬ 
responding number for the Ising energy E is just in the or¬ 
der of thousands. Therefore, the WL simulations for multi¬ 
energy variables cost significantly more in computational time 
than the one-variable case does (see also Refs. 02111511371). 
For instance, our computation of V(E, N ) for the system with 
L = 40 takes about six months on a 3.3 GHz Xeon E3 proces¬ 
sor. While the WL procedures that we consider here is stan¬ 
dard and essentially serial, extending the recently suggested 
broad kernel update method l37l l38l and massively parallel 
algorithm [39j [40| to a multiparameter system may help to 
reduce the issue of long computation time. 

Once r (E,N) is obtained from the WL procedures, it is 
straightforward to calculate the canonical average of a ther¬ 
modynamic observable O = 0(E. N ) at given T and A as 

{O) = \ Y, °( E ’ N ) T ( E ’ N > N e MPE)- (10) 

E,N 

Similarly, one can also define the moment of microcanonical 



FIG. 1. (Color online) Phase diagram of the two-dimensional spin-1 
Blume-Capel model around the tricritical point in the plane of tem¬ 
perature T and crystal field A. The transition line and tricritical point 
are determined from the mixing-field analysis with the calculations 
using the Wang-Landau density of states. The phase-transition line 
plotted here is determined in the limit of infinite L from the extrap¬ 
olation of the pseudotransition points obtained for the systems with 
different sizes up to L = 48 (for example, see Fig. f3}. The transition 
points previously available from the literature I19N331 I42H are given 
for comparison. The statistical error was unspecified for the data in 
Ref. [33], and the error bars of the data in Ref. [19] were given in 
temperature. All our transition points are listed in Table |TT| 

magnetization as 

(M fe ) = 4 '52[(\ m \)E,N} kr ( E ’ N ) zN ex PU3 E )i (H) 

E,N 

where the microcanonical magnetization (|to|)b,jv is an aver¬ 
age of \m\ = L~ d | Si| for a given (E, N) which can be 
measured simultaneously with the WL sampling 11911441 . In 
practice, the microcanonical average is performed in the last 
stage of the iterations with the smallest / where the density of 
states is saturated. In our simulations, the random walk done 
for convergence of the microcanonical magnetization is typ¬ 
ically in the order of a thousand flatness inquiries; however, 
we find that the estimation of (\iti\)e,n is still numerically af¬ 
fordable for L < 40 within the limited computational time. 
We calculate the susceptibility and fourth-order cumulant of 
microcanonical magnetization by using Eq. {TT}. While the 
moment of microcanonical magnetization may quantitatively 
differ from the genuine canonical counterpart, our finite-size- 
scaling analysis given in the later sections shows that it is still 
very useful for the estimate of the first-order-transition points 
and, more importantly, it shares the same universal behavior 
anticipated at the tricritical point. 

m. PHASE DIAGRAM OF THE BLUME-CAPEL MODEL 

In this section, we particularly focus on the area of the 
phase diagram of the BC model at large crystal fields just 
below A = 2. It is known that the first-order transitions 
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FIG. 2. (Color online) Specific-heat anomaly appearing with the 
first-order transition, (a) Double-peak structure observed in the spe¬ 
cific heat at A = 1.994. While the diverging peak at the lower tem¬ 
perature is associated with the phase transition, the anomalous broad 
one at the higher temperature does not scale with the system size. 
The inset indicates that the smaller systems cannot detect the correct 
structure, (b) Corresponding density of nonzero spins (n) and en¬ 
tropy per site s(T). The zero-spin state is dominant right after the 
transition and then starts to get thermally excited to the nonzero-spin 
states distributed with increasing entropy. 


FIG. 3. Finding the first-order transition points by using the method 
of field mixing, (a) Symmetric double peaks in the distribution of 
the field-conjugate variable Q at the phase coexistence A = A* L 
established at temperature T = 0.5. The distribution is shifted and 
rescaled for visualization of systems with different sizes, (b) Scaling 
behavior of AJ, with system size L. The extrapolation in the limit 
of L —^ oo determines the transition point A^. The fourth-order 
cumulant of microcanonical magnetization [7™ in panel (c) also find 
its crossing point very close to the estimated transition point. 


dominates in this area, however, the detailed plot of the first- 
order transition line is not available yet. In Fig. |T| we present 
the transition line that we obtain as a function of temperature 
and crystal field, where the estimated location of the tricritical 
point is also specified. While our estimation of the tricritical 
point is in good agreement with the previous numerical re¬ 
sults G2 ESI E31123 B3. the transition points available from 
the literature fl9l (33 1 shows a significant deviation from the 
first-order transition line that we identify with the WL method. 
The finite-size scaling for the extrapolation is performed with 
system sizes L < 48, and these large-scale calculations are 
also crucial to reveal the valid physics of the specific heat oc¬ 
curring deep in the first-order area. 

Our Wang-Landau calculations also reveal an anomalous 
double-peak structure of the specific heat at large crystal field, 
yet in the first-order-transition area. Figure [2] displays the 
structure of the specific heat at A = 1.994 where the broad 
anomaly emerges above the sharp divergence of the first-order 
transition. This anomaly does not scale with system size and 
is associated with the Schottky-like mechanism. We find that, 
at the first-order transition, the population of nonzero spins 
sharply drops on the disordered side. Since the zero spins 
dominate at this stage, similar to the Schottky anomaly, the 
energy barrier for the excitation of the nonzero-spin states is 
mainly from the crystal field A. While this anomaly could be 


anticipated for very strong anisotropy with A / J 1, it was 

not known whether it could appear together with the phase 
transition. The previous WL calculations for the specific heat 
are limited to L < 16 |jT9l which we find cannot detect the 
double-peak structure [see the inset of Fig.[2ja)]. In our cal¬ 
culation, for system sizes up to L = 48, the double-peak struc¬ 
ture of the specific heat is visible for A > 1.99. 


A. First-order phase transitions 

We determine the first-order phase transition line from the 
symmetry condition for phase coexistence in the probability 
distribution function Pl(Q)- The energy variable Q is con¬ 
jugate to the scaling field A across the phase transition, and 
thus one may expect a symmetric doubly peaked Pl{Q) at 
the transition point in analogy with the probability distribu¬ 
tion of the order parameter in conventional first-order transi¬ 
tions (34). For our systems with finite sizes, we search for a 
size-dependent pseudotransition point and mixing parameter 
at which the symmetric doubly peaked Pl(Q) emerges. 

Figure |3ja) illustrates the symmetric probability distribu¬ 
tion with the double peaks found for the phase coexistence at 
a given temperature. In the graphical search for the symmetry 
to find the pseudotransition point A j and mixing parameter s, 
a practical difficulty lies in discriminating the shape of the dis- 
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FIG. 4. Probability distribution function of the field-conjugate vari¬ 
able Q along the transition line. The calculations for the largest avail¬ 
able system with L = 48 are shown. 

tribution, which actually is the histogram of the discrete values 
of Q constructed with finite bin size. For the optimized iden¬ 
tification of the symmetry, we compare various local statistics 
of the double peaks, including populations and heights. 

The search occurs in three main steps. First, for a given T, 
we graphically search for a set of A and s that roughly gives 
double peaks in Pl{Q) where Q is normalized for Q to have 
zero average and unit variance. Second, starting from these 
initial values, A is fine-tuned to meet the equal population 
condition by measuring the difference in population below 
and above Q = 0. In this step, we also check the symmetry of 
the local averages measured for the parts below and above the 
zero point, which we find comes along with the equal popula¬ 
tion condition. Note that this step is independent of any graph¬ 
ical visualization and therefore allows the high-resolution de¬ 
termination of Aj in the WL approach. Then, with A being 
fixed, the mixing parameter s can be determined graphically 
for the condition of equal height of the double peaks. We find 
that the two peaks are well separated in the first-order transi¬ 
tion region as explicitly shown in Fig. [4] and the tuning of s 
mainly changes the peak heights without disturbing the equal 
population condition. In our numerical implementation, we 
determine the pseudotransition point A j when the difference 
in population and height is minimized within the search step 
of 1CP 6 in A for T < 0.6 and 10 -5 otherwise to get enough 
resolution for size scaling. 

We obtain the transition point A^ from the extrapolation 
of the pseudotransition point Aj in the limit of L —>■ oc . For 
the area of the phase diagram shown in Fig. [I] we find that A j 
shows the scaling behavior 

A* L = A* 00 + bL~ 2 , (12) 

where b is a fitting parameter, which agrees well with the L~ d 
scaling generally expected for the first-order phase transitions. 
The scaling behavior around T = 0.5 is presented in Fig.[3]b), 
for example. We also examine the crossing point of the 
fourth-order cumulant of the microcanonical magnetization 
[7™ = 1 — (|to| 4 )/3(|to| 2 ) 2 measured for systems with dif¬ 
ferent sizes. We find that the crossing of E7™ is in good agree¬ 
ment with the transition point A^ obtained from the analysis 


T 

A* 

^OO 

crossing of [/£* 

order of transition 

0.3 

1.99960(1) 

1.99960 

first 

0.32 

1.99933(1) 

1.99932 

first 

0.34 

1.99895(1) 

1.99894 

first 

0.36 

1.99842(1) 

1.99842 

first 

0.38 

1.99772(1) 

1.99772 

first 

0.40 

1.99681(1) 

1.99681 

first 

0.42 

1.99566(1) 

1.99566 

first 

0.44 

1.99423(1) 

1.99423 

first 

0.46 

1.99248(1) 

1.99248 

first 

0.48 

1.99038(1) 

1.99038 

first 

0.5 

1.98789(1) 

1.98789 

first 

0.52 

1.98496(1) 

1.98496 

first 

0.54 

1.98157(1) 

1.98157 

first 

0.56 

1.97766(1) 

1.97766 

first 

0.58 

1.97323(1) 

1.97322 

first 

0.59 

1.97080(1) 

1.97077 

first 

0.595 

1.96953(1) 

1.96949 

first 

0.6 

1.96825(1) 

1.96817 

first 

0.605 

1.96690(1) 

1.96681 

first 

0.608 

1.96604(1) 

1.96597 

tricritical point 

0.61 

1.96550(1) 

1.96541 

second 

0.615 

1.96412(1) 

1.96399 

second 

0.62 

1.96270(1) 

1.96253 

second 

0.625 

1.96125(2) 

1.96106 

second 

0.63 

1.95980(5) 

1.95954 

second 

0.64 

1.9565(1) 

1.95647 

second 

0.65 

1.9534(1) 

1.95331 

second 

0.66 

1.9501(1) 

1.95006 

second 


TABLE II. Estimated transition points. The extrapolated values of 
A^, obtained from the size scaling in Eq. ( |12| > and the crossing points 
of the fourth-order cumulant U™ are given. 


of the probability distribution [for instance, see Fig.[3jc)]. The 
difference between the two different approaches is observed to 
be about 10 -4 for 0.58 < T < 0.64 beyond the estimated er¬ 
rors which could be further improved by averaging over many 
samples of the WL density of states. The estimated transition 
points are listed in Table. [II] 

The scaling behavior in Eq. ( [T2| certainly supports the first- 
order characteristics of the transition occurring in the area of 
the phase diagram that we are after. Within our data obtained 
for systems with sizes up to L = 48, we have not found 
a quantifiable change of the scaling behavior which, on the 
other hand, one may expect to see above the tricritical point 
of the BC model where the second-order transition should 
emerge. However, in the probability distribution shown in 
Fig.Uj we find that the positions of the double peaks in /\ ( Q) 
get closer as the temperature increases. Above T = 0.64, the 
peaks start to merge together in the larger systems, which im¬ 
plies that the character of the transition indeed alters. 
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B. Tricritical Point 

We determine the precise location of the tricritical point 
from the scale-invariant universal form of the probability dis¬ 
tribution function Pl(Q), as indicated in Eq. (|8j. The scale 
invariance at the tricritical point can be conveniently indicated 
by a size-independent crossing point of the fourth-order cu- 
mulant 


U% = 1 - 


(Q 4 ) 

3(Q 2 ) 2 


(13) 


for the field-conjugate variable Q normalized to have zero av¬ 
erage and unit variance mmm. 

We identify the tricritical temperature as T tc = 0.6080(1) 
from the location of the crossing point of the fourth-order cu- 
mulant t/® along the transition line as shown in Fig. [ 5 J a). 
Note that the transition line here is for finite L: namely, the 
line of the pseudotransition points A * L {T) that we have de¬ 
termined for the phase coexistence. The error estimation is 
only graphical since the our calculation is based on a sin¬ 
gle sample of the WL density of states. We estimate the 
tricritical crystal field as A tc = 1.9660(1) from the extrap¬ 
olation of the pseudotransition point A) and also from the 
crossing point of the fourth-order cumulant (7™ measured at 
T = T tc [ see Figs. [5fb) and [5jc)] . Our estimate of the 
tricritical point, T tc = 0.6080(1) and A tc = 1.9660(1), 
is in very good agreement with the previous results for the 
spin-1 Blume-Capel model in square lattices, which provide 
T tc = 0.610(5) and A tc = 1.965(5) (33), T tc = 0.608(1) 
and A tc = 1.9665(3) (34), T tc = 0.609(4) and A tc = 
1.965(5) ED, T tc = 0.609(3) and A tc = 1.966(2) OH, and 
very recently T tc = 0.608(1) and A tc = 1.9665(3) l36l . 


IV. TRICRITICAL SCALING BEHAVIOR 


In this section, we present the three different forms of 
finite-size-scaling analysis that we perform to determine the 
tricritical eigenvalue exponents. The thermal exponent y t 
is extracted from the probability distribution function of the 
field-conjugate variable Q at the tricritical point. The scal¬ 
ing of the fourth-order cumulant [ 7 ® along the transition 
line is examined to obtain the next-to-leading thermal expo¬ 
nent y g . Finally, we perform the phenomenological finite- 
size-scaling analysis with thermodynamic quantities includ¬ 
ing specific heat, compressibility, susceptibility, magnetiza¬ 
tion to measure the thermal and magnetic exponents y t and 
Vh- 


A. Distribution of the field-conjugate variable 

We examine the tricritical thermal exponent y t from the 
probability distribution function given in Eq. <§• Precisely 
at the tricritical point, T = T tc and A = A * L {T tc ), the distri¬ 
bution function for the relevant field-conjugate variable Q can 




FIG. 5. Location of the tricritical point, (a) The tricritical tempera¬ 
ture T tc ~ 0.6080 is determined at the crossing point of the fourth- 
order cumulant of the field-conjugate variable [7® along the transi¬ 
tion line A = A/(T). The extrapolation of the transition points 
Ai, (Ttc) in (b) and the crossing point of the fourth-order cumulant 
of microcanonical magnetization 17™ in (c) provide the estimation of 
the tricritical crystal field as A tc ~ 1.9660(1). 


be reduced into the simple finite-size-scaling ansatz (34) as 

Pl(Q) = L d -ytp* Q (L d ~y'Q), (14) 

where pg is a universal function and the dimension is given 
as d = 2 for square lattices. 

Figure [6| a) presents our finite-size-scaling analysis for the 
probability distribution with the tricritical thermal exponent 
yt = 1.80(1), showing the data of Pl(Q) falling well onto 
a single curve. In particular, the lines for L = 40 and 48 
can hardly be distinguished in the plot because of the almost 
perfect overlap. The possible error in this estimate with the 
shape of the distribution mainly originates from the discrete 
nature of Q which affects the visualization of its histogram, 
particularly in small systems, and thus can cause ambiguity in 
the graphical determination of the mixing parameter. 

Our estimation of y t — 1.80 numerically confirms the ex¬ 
act conjecture y t = 9/5 within the graphical identification. 
The data collapse of Pl(Q) for systems with different sizes 
L = 16 to 48 shows good agreement with the previous finite- 
size scaling for the spin fluid model, which was also compared 
for universality with the BC model with size L = 40 041 . In 
principle, one can also attempt to extract the next-to-leading 
exponent y g from the similar finite-size scaling of the proba¬ 
bility distribution function Pl{£) as implied in Eq. (|8|. How¬ 
ever, we find that Pl{£) does not give any meaningful esti¬ 
mate of y g because the distribution is too close to the Gaus¬ 
sian normal distribution, regardless of the system size L. The 
same issue was also reported by the previous work l34l where 
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FIG. 6. (Color online) Finite-size-scaling tests of the field-conjugate 
variable Q for the tricritical thermal exponents, (a) Scaling plots 
of the probability density distribution Pl(Q) with the exponent 
y t = 1.8 at the tricritical temperature T tc = 0.608. The large sys¬ 
tems with L > 40 provide smooth curves (solid lines) falling on a 
universal distribution which also fits well with the data points for the 
smaller systems (symbols), (b) Finite-size scaling of the fourth-order 
cumulant Uf along the transition line. The data points fall onto each 
other very well, given the next-to-leading exponent y g = 0.8. 


y g = 1.03(7) was estimated from the finite-size-scaling test 
of P l (£). 


B. Fourth-order cumulant along the transition line 

Instead, we utilize the fourth-order cumulant ( 7 ® for the 
estimate of the next-to-leading thermal exponent y g . From Eq. 
(J7|! as well as from the scaling hypothesis of the persistence 
length j33l . one can find that the finite-size scaling of Ilf 
along the transition line may follow the scaling form Uf = 
u*[L Vg g\, where u* is a universal function, and the scaling 
field g is the deviation from the tricritical point in the direction 
tangent to the coexistence curve. Moreover, in our observation 
of the data for the phase diagram, it turns out that (/x — nt) is 
almost linearly proportional to (j3 — j3 1 ) along the transition 
line near the tricritical point, which leads to g ~ (T — T tc ) in 
Eq. [4] Therefore, for the explicit finite-size scaling tests, one 
can further simplify the scaling ansatz of Uf as 

E#|a=A£(t) « ti*[i v «(T-T tc )], (15) 


Figure[6|h) shows that our data points of Uf along the tran¬ 
sition line fall perfectly onto the same curve in the test with 
y g = 0.8 for Eq. ( p~5j ). Within the graphical uncertainty, we 
determine the next-to-leading thermal exponent y g = 0.80(1), 
which confirms the exact conjecture y g = 4/5. Our finite- 
size-scaling analysis of (jf can be compared with the finite- 
size-scaling test of the persistence length, which indicates 
y g = 0.80(1) li33l and the estimate made by using the slope of 
the fourth-order cumulant which provides y g = 0.83(5) ll34ll . 


C. Phenomenological finite-size scaling 


In this section, we present the phenomenological finite- 
size-scaling analysis of thermodynamic quantities to deter¬ 
mine the thermal and magnetic exponents y t and ;<//, . This 
approach does not directly rely on the field-conjugate variable 
Q and its probability distribution function. Therefore, it is free 
from the explicit dependence of the mixing parameter and the 
histogram-visualization issue for the discrete data of Q. 

We consider susceptibility, magnetization, specific heat, 
and compressibility as the thermodynamic quantities to be ex¬ 
amined for our finite-size-scaling analysis. The susceptibility 
X = (L d /T)((\m\ 2 ) — (|m|) 2 ) and the magnetization (|m|) 
are estimated with the microcanonical magnetization by using 
Eq. (111. The specific heat c = (L d /T 2 )((e 2 ) — (e) 2 ) and 
the compressibility kt = (L d /T)((n 2 ) — (n) 2 )/(n) 2 are re¬ 
lated to the fluctuations of the energy E and the number N 
of nonzero spins. With the WL density of states being sam¬ 
pled with high accuracy, one can freely access these thermo¬ 
dynamic variables at any temperature and crystal field. 

Figures [7] and [8] show our finite-size-scaling analysis of the 
thermodynamic quantities for two different choices of an ap¬ 
propriate scaling axis. First, we choose to perform the finite- 
size scaling along the fugacity axis selected from the natural 
variables of the grand partition function. With the tempera¬ 
ture fixed at T = T tc , the scaling variable can be expressed as 
x = A — Ate■ In this case, the relevant thermodynamic quan¬ 
tities are the number fluctuations, susceptibility, and magne¬ 
tization, while the specific heat is discarded for our choice 
of the scaling test with fixed T. The corresponding scaling 
ansatz can be written as 


™) 2 K T = L“‘ / " i r(iL 1/i '*), 

x = L ^ /Vt x°( xL 1/Vt ), 

(|m|) = L~M Vt M°(xL 1 / Vt ), 


(16) 

(17) 

(18) 


where A f°, x°> and Ad " are universal functions. In compar¬ 
ison with Eq. ([7J, one can also obtain the relations between 
the conventional exponents, v t , at, fit, and 7 t , through the 
tricritical eigenvalue exponents y t and ///, as 


a-t/vt = ~d + 2y t . 
~Pt/vt = -d + yh, 
It I ft = -d + 2y h . 


(19) 

( 20 ) 
( 21 ) 


where the constraint A = A * L {T) ensures that it is along the Provided the hyperscaling identity v t d = 2 — a t , the thermal 
transition line for a system with finite size L. exponents are simply related as y t = 1 /Vf. 
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-1.5 -1 -0.5 0 0.5 1 1.5 

(A-A tc ) l3 

FIG. 7. (Color online) Tricritical behavior along the crystal field 
axis at the tricritical temperature. The finite-size-scaling analysis of 
(a) number fluctuations {ti) 2 kt, (b) susceptibility and (c) mag¬ 
netization (|m|) is performed to determine the tricritical exponents. 
While the WL method guarantees high enough resolution to plot the 
data as continuous curves, the points in low resolution (L < 24) are 
also given for visualization of finite-size scaling. The ratios at /vt 
and 7 */vt are determined from the power-law fits of the maxima of 
(n) 2 Kr and respectively. Each scaling plot with the estimated 
exponent shows the excellent collapse of the data points falling onto 
a single curve. The tricritical eigenvalue exponents are deduced as 
yt = 1.804 and yh = 1.925 from at/vt and 7 t./vt- 


The thermal exponent y t can be easily extracted from the 
maxima of ( n) 2 KT which scales as ( n) 2 KT oc L at / Vt . Fig¬ 
ure [7ja) shows the power-law fit of the maxima, providing 
the estimate of 07/17 = 1.608. This ratio of the exponents 
can be directly converted into the tricritical thermal exponent 
as y t = 1 /v t = 1.804 which turns out to be very close to 
the exact conjecture y t = 9/5. The full finite-size-scaling 
ansatz for (ri) 2 KT is also examined with the estimated expo¬ 
nents 07/77 = 1.608 and \/v t = 1.804, showing the excellent 
collapse of the data curves falling onto a single line, as shown 
in Fig. [7} a). 

We estimate the magnetic exponent ;<//, through the similar 
analysis for the susceptibility of which maxima scales as \ cx 
U <l l Vl . From the power-law fit shown in Fig. 0 b>. we find 
out 7 t/vt = 1-850, and this ratio is directly converted into 
the tricritical magnetic exponent y^ = 1.925 which precisely 
agrees with the exact conjecture yh = 77/40. Figure [7jb) 



-1 -0.5 0 0.5 1 

(T-T tc ) L 1/Vt 

FIG. 8 . (Color online) Tricritical behavior along the temperature 
axis. The fugacity is fixed at Ins = A/T = Atc/Tt c . The finite- 
size scaling analysis of (a) specific heat c, (b) susceptibility 7 , and 
(c) magnetization m is presented. The tricritical exponents are deter¬ 
mined by the same procedures used in Fig. [7] The resulting scaling 
plots show the excellent collapse of the data points falling onto a sin¬ 
gle curve, where the corresponding tricritical eigenvalue exponents 
are deduced to be yt = 1.809 and yh = 1.9275 from at/vt and 

'It/vt. 


shows the data perfectly falling onto a single curve in the test 
of the finite-size-scaling ansatz, confirming the accuracy of 
our estimate of the magnetic exponent. For the magnetization, 
while fit/vt can be directly determined by the obtained 7 t /v t 
by using the scaling relations through yh, we also examine 
the finite-size scaling ansatz of (to) for explicit confirmation, 
where we find the excellent collapse of the data curves falling 
onto a single line, as indicated in Fig.[7Jc) 

One the other hand, we perform another estimate of the 
tricritical exponents by choosing the T axis for the similar 
finite-size-scaling analysis. The fugacity z is now fixed at 
biz = A/T = A t c /T tc , and thus the scaling variable is 
given as x = T — T tc . In this case, the relevant thermo¬ 
dynamic quantity for finite-size scaling includes the specific 
heat; namely, the energy fluctuations, instead of the number 
fluctuations. Although, the finite-size-scaling ansatz for the 
specific heat c can be written similarly as 

c = L at l Vt C°(xl}l Vt ), 


( 22 ) 
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where C° is a universal function. The same scaling relation 
between at /vt and yt holds for the specific heat as well. By 
applying the same procedures as done for the earlier finite-size 
scaling in the A-axis, here we estimate the tricritical expo¬ 
nents as y t = 1.809 and yh = 1.9275 on the T axis, as shown 
in Fig. [8] While the estimate of the tricritical exponents on 
the T axis are slightly different from those estimated in the 
finite-size scaling on the A axis, both estimations are still in 
very good agreement with the exact conjectures, y t = 9/5 and 
y t = 77/40. The source of the discrepancy found between 
the two estimates may originate from the possibility that the 
error in locating the tricritical point propagates differently in 
our two choices of the scaling and fixed variables in the phe¬ 
nomenological finite-size-scaling analysis. 

Finally, from the different forms of finite-size scaling that 
we have performed so far in this section, we can write the 
tricritical eigenvalue exponents of the BC model as 

yt = 1.804(5), y g = 0.80(1), y h = 1.925(3), 

showing very good agreement with the exact conjectures, 
yt = 9/5, y g = 4/5, and yh = 77/40, respectively. The esti¬ 
mated errors are mainly from the slight difference between the 
values observed in the different approaches of finite-size scal¬ 
ing. The comparison with the previous works using different 
numerical methods are also listed in l ahle Q] 


V. CONCLUSIONS 

In conclusions, we have demonstrated the effectiveness of 
the Wang-Landau method in the finite-size-scaling analysis 
for tricritical behavior within the spin-1 Blume-Capel model 
in two dimensions. The significance of our results is two-fold. 
First, we have constructed the detailed line of first-order tran¬ 
sitions, completing the previously-less-explored area of the 
phase diagram at low temperatures, which is hardly accessi¬ 


ble in conventional Monte Carlo simulations. In the area of 
large crystal fields very close to A = 2, we have found a 
double-peak structure in the specific heat where the Schottky- 
like anomaly is observed above the first-order-transition tem¬ 
perature. Second, through the various forms of the finite- 
size-scaling analysis, we have successfully estimated the tri¬ 
critical point as T tc ~ 0.6080 and A tc — 1.9660 and the 
tricritical exponents as y t = 1.804(5), y g = 0.80(1), and 
yh = 1.925(3). In particular, our high-resolution analysis of 
the phenomenological finite-size scaling takes a great advan¬ 
tage from the Wang-Landau methods, granting unrestricted 
access to the values of temperatures and crystal fields. 

The performance of the Wang-Landau method may depend 
on its practical limit in the system size which is still much 
smaller than those accessible in conventional methods. The 
large computational resource requirement is indeed one of 
the biggest obstacles that the Wang-Landau method should 
overcome to show its effectiveness in challenging problems 
of phase transitions. We have shown that, within the limit of 
our computational resources, the standard Wang-Landau al¬ 
gorithm now allows us to simulate the Blume-Capel model 
with sizes up to 48 x 48 sites, which provide excellent finite- 
size scaling for the tricritical behavior. Our demonstration 
suggests that, with increasing computational power and po¬ 
tential support from more advanced techniques such as the 
recently suggested parallel algorithm for scalability Ii39l l40ii. 
the Wang-Landau method may provide a promising tool of 
high-precision numerics for multicritical phenomena. 
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